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1 .  Introduction 


A  common  application  of  FASCODE l>2  calculations  is  a  comparison  to  measured 
spectra.  In  order  to  make  a  proper  comparison,  the  calculated  spectrum  must  be  convolved,  or 
smoothed,  with  the  spectral  response  function  appropriate  to  the  instrument  being  modeled. 
Frequently  the  measurements  are  made  with  a  Fourier  Transform  Spectrometer  (FTS)  for  which 
the  spectral  response  function  is  known  and  which  may  be  modified  during  data  reduction1 2 3.  By 
apodizing  the  interferogram,  the  shape  and  resolution  of  the  spectral  response  function  can  be 
controlled. 

In  previous  versions  of  FASCODE,  spectral  smoothing  or  "scanning"  was  performed  by 
an  actual  convolution  in  the  spectral  domain.  The  scanning  function  is  calculated  over  some 
finite  spectral  extent—typically  some  number  of  halfwidths  or  zeros  from  the  line  center— and  is 
convolved  with  the  calculated  spectrum.  The  accuracy  of  the  result  depends  strongly  on  the 
spectral  extent:  for  a  function  like  sine,  whose  sidelobes  fall  off  slowly,  the  scanning  function 
must  be  calculated  out  to  more  than  100  halfwidths  to  achieve  accuracies  at  the  0.2  percent  level. 
However,  the  computational  time  varies  linearly  with  the  spectral  extent  of  the  scanning  function 
so  that  accurate  calculations  with  functions  like  sine  can  be  expensive. 

We  have  added  the  capability  to  FASCODE  to  model  the  spectral  response  function 
associated  with  an  FTS  by  using  Fourier  transforms.  This  technique  directly  mimics  the 
operation  of  an  FTS;  the  calculated  spectrum  is  transformed  into  an  "interferogram ",  "apodized". 
and  then  transformed  back  to  the  spectral  domain.  This  technique  is  more  accurate  than  the 
convolution  technique  and  in  many  cases,  more  efficient.  The  user  is  given  the  choice  of  a 
number  of  commonly  used  apodization  functions  to  define  tne  shape  of  the  scanning  function, 
and  may  specify  the  resolution  in  terms  of  either  the  half  width  of  the  scanning  function  or  the 
optical  path  difference  of  an  equivalent  interferometer. 


1  Clough,  S.,  A.,  F.  X.  Kneizys,  L.  S.  Rothman,  and  W.  O.  Gallery,  1981, Atmospheric  spectral 

transmittance  and  radiance:  FASCOD1B,  proc.  of  Soc.  Photo.  Opt.  Instrm.  Eng..  277. 
152-166. 

2  Clough,  S.  A.,  F.  X.  Kneizys,  G.  P.  Anderson,  E.  P.  Shettle,  J.  H.  Chetwynd.  L.  W.  Abreu.  and 

L.  A.  Hall,  1989,  FASCOD3:  Spectral  Simulation,  in  IRS  '88:  Current  Problems  in 
Atmospheric  Radiation.  J.  Lenoble  and  J.  F.  Gelyn  (Eds.),  A.  Deepak  Pub. 

3  Bell,  R.  J.,  1972,  Introductory  Fourier  Transform  Spectroscopy.  Academic  Press,  New'  York 
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2.  Theory 

We  will  begin  by  briefly  describing  the  mathematics  of  Fourier  Transform  Spectroscopy 
and  then  make  the  connection  with  the  present  work. 

2.1  Fourier  Transform  Spectroscopy 

The  quantity  measured  by  an  interferometer  is  the  interferogram  I(x)  as  a  function  of  the 
optical  path  difference  x.  For  an  ideal  spectrometer,  /  is  related  to  the  incident  spectrum  S(v) 
through  the  Fourier  transform  J: 


+oo 


I(x)  =  f(S)  =  f  5(o)  exp(-2jtirv))  do 

-oc 


(1) 


In  practice,  I(x)  is  only  measured  out  to  L,  the  maximum  optical  path  difference  of  the 
interferometer.  The  spectrum  S'(o)  recovered  from  the  interferogram  is  obtained  from: 

+L 

5'(u)  =  ‘/(AT)  =  fA(x)  l(x)  exp(-2mi))  d.r  (2) 

■L 

where  A(x)  is  the  apodization  function  applied  to  the  interferogram  to  control  the  shape  of  the 
scanning  function. 

The  convolution  of  two  functions  5  and  R  is  defined  by  the  following  expression; 

SXo)  =  R  ☆  5  =  J/?(o)  5(o’-u)  dt)'  (3) 


where  the  symbol  ☆  represents  convolution.  A  fundamental  theorem  of  Fourier  transforms 
states  that  the  convolution  of  two  functions  equals  the  transform  of  the  product  of  the  transforms 
of  the  individual  functions.  Applying  this  theorem  to  Eq.  2  gives: 

S'(u)  =  ?(A)  ☆  W)  =R(v)  ☆  5(u)  (4) 

where  R(v)  =  /(A(x)).  R(v)  is  the  scanning  function  associated  with  the  apodization  function 
A(x). 

If  no  apodization  is  applied,  then  A{x)  is  effectively  a  rectangle  of  width  equal  twice  the 
optical  path  difference.  The  scanning  function  associated  with  this  apodization  is  the  sine 
function  defined  as: 


sinc(o)  =  sinf2nioL)/(2KioL) 


(5) 


This  function  is  characterized  by  relatively  large  side  lobes  which  fall  off  slowly,  ft  is  possible 
to  trade  off  resolution  against  smaller,  more  rapidly  decaying  side  lobes  by  choosing  different 
apodization  functions.  Figure  1  shows  five  common  apodization/scanning  function  pairs.  These 
figures  correspond  to  the  HIS4  interferometer  with  a  maximum  path  length  difference  of 
1.3735  cm.  This  instrument  was  chosen  as  a  representative  Fourier  transform  spectrometer. 

In  practice,  interferograms  are  discretely  sampled  on  a  grid  xx  =  iAt,  i=  0  to  N-l,  where 
Ax  is  the  sampling  interval  and  N  is  the  number  of  points.  (Since  the  ideal  interferogram  is 
symmetric  around  zero,  only  the  positive  half  need  be  retained.)  The  continuous  Fourier 
transform  in  Eq.  2  is  replaced  with  the  discrete  transform: 

N-l 

S’(vj)  =  ‘I (AD  =  Z'3(x()  I(x j)  expl-lTtLVjUj)  (6) 

i=0 

The  frequency  grid  is  given  by  uj  =  j  Au,  j=  0  to  N-l,  where  Au  =  1/(2 L).  The  Nvquest 
frequency,  umax  =  (N-l)Au  =  l/(2Av).  is  the  highest  frequency  that  can  be  properly  sampled.  If 
frequencies  higher  than  this  are  present  in  the  signal,  then  they  will  be  aliased  down  to  lower 
frequencies  in  the  recovered  spectrum  and  the  spectrum  will  be  distorted. 

2.2  Spectral  Smoothing 

In  order  to  compare  a  FASCODE  spectrum  to  a  measured  spectrum,  the  calculated 
spectrum  must  be  convolved  with  the  scanning  function  appropriate  to  the  instrument  being 
modeled.  This  convolution  has  been  implemented  using  Fourier  transforms: 

S(v)  =  ‘f(  J(R)  •  -KS)  )  =  R  ☆  5  (  7 1 

where  S  is  the  monochromatic  spectrum,  R  is  the  scanning  function,  and  S'  is  the  smoothed 
spectrum.  By  analogy  to  Fourier  transform  spectroscopy,  ‘J( S)  is  the  interferogram  and  f(R)  is 
the  apodization  function. 

Note,  however,  that  the  calculated  spectrum  is  limited  from  U|  to  vh-  Before 
transforming  the  spectrum,  it  is  shifted  in  frequency  down  to  the  range  0  to  x>2-V\.  In  the 
interferogram  domain,  we  now  have  Ax  =  l/(2(u;>-ui))  and  L  =  l/(2Ao). 

However,  a  significant  difference  exists  between  smoothing  a  calculated  spectrum  and 
apodizing  an  interferogram,  having  to  do  with  edge  effects.  As  a  consequence  of  using  discrete 
Fourier  transforms,  the  measured  spectrum  must  be  thought  of  as  repeating  infinitely  in  the 
positive  and  negative  directions,  with  a  period  of  umax  =  1/Ar.  Near  the  edges  of  the  spectrum, 
at  0  and  umax,  the  wings  of  the  scanning  function  encounter  the  repeated  spectra. 


4  Revercomb,  H.  E.,  H.  Buijs,  H.  B.  Howell.  D.  D.  LaPortc,  W.  L.  Smith,  and  L.  A.  Sromovsky. 
1988,  "Radiometric  Calibration  of  IR  Fourier  Transform  Spectrometers:  Solution  to  a 
Problem  with  the  High-Resolution  Interferometer  Sounder,"  Appl.  Opt,  27.  p32IO 
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Figure  1.  Five  Common  Scanning  Function/Apodization  Function  Pairs:  Sinc/Rectangle. 
SincVTriangle,  Beer/(l-(.r/C)2)2,  Hamming,  and  Hanning.  The  maximum  displacement  of 
1.3735  cm  corresponds  to  the  HIS  instrument. 
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The  calculated  spectrum  is  limited  from  \>(  to  ut.  or.  after  shifting,  from  0  to  iq-in.  This 
spectrum  repeats,  with  a  period  of  so  that  effectively  wraps  around  to  Uj.  Near  the 

edges  of  this  spectrum,  the  wings  of  the  scanning  function  encounter  a  spectrum  different  from 
that  seen  in  the  measured  spectrum.  This  effect  is  especially  severe  when  modeling  unapodized 
spectra  with  the  sine  scanning  function,  whose  side  lobes  are  particularly  large.  In  order  to 
minimize  this  problem,  the  range  of  the  calculated  spectrum  (Uj  to  vhi  should  be  made 
significantly  larger  than  the  actual  range  of  interest. 

The  seven  scanning  functions  included  in  FFTSCAN  are  listed  in  Table  1  along  with  the 
corresponding  apodization  function.  For  an  FTS,  the  most  convenient  measure  of  the  width  of 
the  scanning  function  is  in  terms  of  the  parameter  a  which  equals  the  reciprocal  of  the  maximum 
optical  path  difference  L  of  an  equivalent  interferometer,  a  is  approximately  the  resolution  of 
the  apodized  spectra,  according  to  the  Rayleigh  criterion.  It  is  also  common  to  characterize  the 
resolution  of  an  instrument  in  terms  of  the  halfwidth  at  half  maximum  (HVVHM)  and  the 
distance  to  the  first  zero  (FZ).  The  ratios  u/HWHM  and  n/FZ  are  listed  in  Table  1  and  can  be 
used  to  convert  from  one  form  to  another. 

Theoretically,  the  "apodization"  function  ‘f\R)  can  bo  computed  either  us  the  Fourier 
transform  of  the  scanning  function  or  analytically  from  the  function  listed  in  Tabic  I  The 
analytic  method  is  more  efficient  but  it  assumes  that  the  scanning  function  is  infinite  in  extent. 
Because  discrete  transforms  are  used,  the  extent  of  the  scanning  function  is  limited  to  u->  -  tq. 
The  discrete  transform  of  the  scanning  function  over  this  extent  will  not  exactly  equal  the 
apodization  function,  and  if  the  halfwidth  of  the  scanning  function  approaches  the  range 
-  Ui,  then  the  difference  can  become  significant.  This  effect  is  illustrated  in  Figure  2.  Here 
the  scanning  function  is  a  sine2  with  a  HWHM  or  0.3227  cm'1,  which  corresponds  to  the 
resolution  of  the  HIS  instrument.  This  scanning  function  is  to  be  applied  to  a  spectrum  with  an 
extent  of  12.9  cm1,  which  corresponds  to  40  hulfwidths  (HWHM/d)^  -  iq)  =  Cr  in  Table  1 ). 
The  analytic  apodization  function  is  a  triangle  of  base  1.376  cm.  The  plot  labeled  "Error  in  Scan 
Function"  shows  the  difference  between  the  analytic  sine2  function  and  the  FFT  of  the  triangular 
apodization  function.  Similarly,  the  plot  labeled  "Error  in  the  Apodization  Function"  shows  the 
difference  between  the  analytic  triangle  function  and  the  FFT  of  the  sine2  scanning  function. 
This  effect  is  greatest  for  functions  with  sharp  edges,  e  g.  the  rectangle  and  triangle.  The 
column  Cr  in  Table  1  lists  for  each  function  the  critical  ratio  of  (\)~>  -  iq  i/HWHM  at  w  hich  the 
maximum  error  in  the  scan  function  due  to  this  effect  becomes  approximately  0.05  percent. 

In  view  of  this  fact,  the  program  has  been  designed  with  the  capability  to  calculate  the 
apodization  function  by  either  method.  By  default,  the  program  chooses  which  method  to  use 
based  upon  the  ratio  (Ot  -  0{)/HWHM.  This  default  can  be  overridden  {see  user  instructions.) 
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Table  1.  Scanning  Function/Apodization  Function  Pairs 


# 

Name 

Stalling  Function 

Apodization  Function 

a/HWHM 

a/FZ 

CR 

I. 

Tnangie/Sir  w- 

1  -  W a,  lul  1  a.  0:  lul  >  a 

(sin(7wu)/(  rutin- 

2.0 

1 .0 

40 

T 

Gauss/  Gauss 

exp<-0.5  cu/«i“) 

exp(-27t  fa..-)-) 

0.840322 

iNAi 

to 

3- 

Sinc~/Triangle 

(sin(jn)/u)/(!ru/i;))3 

1  -  xu.  I.d  <  1  /a.  0:  l.vl  >  I la 

2,257609 

10 

40 

». 

Sinc/Rectangle 

sin0<)/0<) 

I:  Lrl  <  1  /a.  0:  hi  >  I/a 

3.3  i  4800 

2  0 

160 

5. 

Beer 

J(5/2  ,u)/m<5/2) 

(1  -  (vu)3)3 

2.100669 

20 

6. 

Hamming 

sinc(u)  +  C)  (sinc(i<+ir) + 
sinc(u-Jtl) 

(1  +  2c  |  cost  xvu )  >/(  1  +2c  | ) 

2.195676 

1.0 

20 

7. 

Hanning 

sinc(u)  +  .5  (sinc(u+rc)  + 
sinc(n-Ji)) 

( 1  +  cos(JLra))/2 

2.0 

1.0 

20 

Notes: 

a.  o  =  frequency,  in  cm'1,  x  =  optical  path  difference,  in  cm.  u  =  2~'o/a 

b.  a  -  ML.  where  L  is  the  maximum  optical  path  difference  of  an  equivalent 

interferometer,  a  determines  the  resolution,  or  the  width  of  the  scanning 
function. 

c.  HWHM  the  halfwidth  at  half  maximum  of  the  scanning  function. 

d.  FZ  is  the  distance  from  the  center  of  the  scanning  function  to  the  first  zero. 

e.  CR  =  the  critical  value  of  the  ratio  of  the  extent  of  the  spectrum  ih  -  tq  and  the 

HWHM.  When  this  ratio  is  less  than  CR.  the  apodization  function  is  calculated  as 
the  FFT  of  the  scanning  function.  When  the  ratio  is  greater,  it  is  calculated 
analytically.  (See  the  text.) 

f.  J (5/2,m)/m(3/2)  =  ((3 -u*)  sin(u)  -  3uccMu))/tt~\  J(n,«)  is  the  Bessel  function  ot 
order  n 

g.  c,  =0.428752 

By  necessity,  both  the  interferogram  and  the  recovered  spectrum  from  an  FTS  are  real 
functions.  The  transform  of  a  real  function  is  symmetric  and  that  of  a  symmetric  function  is 
real,  so  that  both  the  interferogram  and  the  spectrum  must  be  symmetric.  If  only  a  one-sided 
interferogram  is  measured,  then  it  must  be  reflected  around  the  origin  to  produce  a  two-sided 
interferogram  before  transforming.  However,  the  convolution  in  Eq.  7  applies  to  any  real 
functions,  symmetric  or  not.  Therefore,  the  calculated  spectrum  S  need  not  be  reflected  to 
produce  a  symmetric  function  before  transforming.  (The  scanning  function  R(\))  is  symmetric.) 
The  "interferogram"  ‘J(S)  will  be  a  complex  function  but  the  convolved  spectrum  S'  will  again 
be  real. 
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Figure  2.  Analytic  Versus  Discrete  Fourier  Transform  Calculations  of:  the  Sine-  Scanning 
Function  and  the  Triangular  Apodization  Function.  The  plot  labeled  "Error  in  Scan  Function"  is 
the  difference  between  the  sine-  function  and  the  FFT  of  the  triangle  while  the  plot  labeled 
"Error  in  Apod  Function"  is  the  difference  between  the  triangle  and  the  FFT  of  the 
sine2  function.  See  the  text  for  further  explaination. 

2.3  Prescanning  With  a  Boxcar 

In  cases  where  the  width  of  the  scanning  function  is  large  compared  to  the  frequency 
spacing  of  the  input  spectrum,  significant  computational  savings  can  be  achieved  by  first  con¬ 
volving  the  spectrum  with  a  rectangle  whose  width  is  small  compared  to  that  of  the  scanning 
function.  In  this  procedure,  M  adjacent  points  are  averaged  to  one  output  point  at  the  mean 
frequency  of  the  M  points.  This  form  of  convolution,  referred  to  here  as  prescanning  with  a 
boxcar,  is  very  fast  and  reduces  the  number  of  output  points  by  a  factor  of  M.  The  resulting 
spectrum  is  then  convolved  with  the  desired  scanning  function  using  Fourier  transforms. 
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Mathematically,  the  procedure  relies  on  the  fact  that  convolution  is  associative.  Let  5  be 
the  input  spectrum,  Rx  be  a  scanning  function  of  halfwidth  at  half  maximum  a,  and  R2  be  a 
rectangle  of  half  width  M  Au/2.  Then: 

S"~(Si?R^)  =S£r(/?1  ☆  £,)  (8) 

Let  r  =  a/(M  A\>/2).  If  r»  1,  that  is,  if  the  rectangle  R ^  is  narrow  compared  to  the  scanning 
function  R  ,  then  R2  ft  Rx  ~  Rx  and  5"  ~  S  =  S  ft  R j.  Here,  "narrow''  is  taken  to  mean  that  /  > 
8.  The  error  introduced  by  this  approximation  will  be  shown  later. 

It  is  also  possible  to  partially  compensate  for  the  error  introduced  by  boxcaring.  This 
procedure  is  as  follows.  Using  the  notation  from  Eq.  4: 

S'  =  !F(S  ft  /?[) 


=  •/(  f(S)  ■  nR |)) 

=  n  'i(S)  ■  :f{R2)  ■  'j(Rx)i  '/(r2)) 

=  n  ns* r2)-  nROf  ‘i(R2))  <9) 

In  the  last  form  of  Eq.  9,  dividing  by  7{R2)  removes  the  effect  of  the  convolution  of  5  with  R2. 
Since  R ^  is  a  narrow  rectangle,  J(R2)  is  a  broad  sine  function.  In  practice,  the  apodized 
"interferogram"  is  divided  by  J{R2)  before  it  is  transformed  back  to  the  spectral  domain. 

The  actual  computation  of  S’  is  performed  using  discrete  transforms,  not  integrals,  and 
the  convolution  shown  as  S  tr  R2  resamples  the  spectrum,  taking  only  every'  M'th  point.  It  is  the 
resampling  which  produces  the  computational  efficiency,  since  it  reduces  the  number  of  points 
in  the  subsequent  Fourier  transforms  by  a  factor  of  M.  However,  because  of  resampling,  the 
deconvolution  process  is  not  exact;  therefore,  it  is  not  possible  to  eliminate  completely  the 
smoothing  effects  of  the  boxcar.  The  error  associated  with  boxcaring  will  be  discussed  next. 

Boxcaring  introduces  two  type  of  errors— frequency  displacement  and  convolution  error. 
The  first  type  of  error— frequency  displacement— is  illustrated  in  Figure  3  which  show's  the  effect 
of  boxcaring  a  delta  function.  As  can  be  seen  in  the  figure,  the  position  of  the  delta  function  can 
be  shifted  by  as  much  as  M  Ai>/2,  depending  on  the  alignment  of  the  resampled  grid  relative  to 
the  delta  function.  A  similar  shift  can  occur  in  a  spectral  feature  of  small  but  finite  width,  as  will 
be  seen  later. 

Boxcaring  also  introduces  an  error  due  to  the  convolution  effect  of  the  boxcar,  which  can 
only  be  partially  compensated  for  by  the  deconvolution  process.  This  error  plus  the  frequency 
displacement  error  are  both  illustrated  in  Figure  4.  In  this  example,  the  (synthetic)  input 
spectrum  is  a  single  Lorenz  line  of  HWHM  =  0.08  cm' 1  on  a  grid  of  0.0 1  cm' 1 .  The  scanning 
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Figure  3.  Schematic  Diagram  Showing  the  Effect  of  Boxcaring  A  Delta  Function.  The  input 
spectrum  containing  a  delta  function  is  boxcared  and  resampled  with  a  value  of  M  =  7.  The 
output  spectrum  contains  a  delta  function  shifted  in  frequency  by  MAu/2. 

function  is  a  sine2  of  HWHM  =  0.3  cm'1.  A  value  of  r  -  8  gives  a  value  of  M  = 
2x0.3/(8x0.01)  =  7  and  a  boxcar  of  half  width  =  7x0.1/2  =  .035  cm'1.  The  output  spectra 
are  therefore  on  a  grid  of  0.07  cm'l.  The  convolved  spectrum  has  been  calculated  without 
boxcaring  and  with  boxcaring  for  4  different  alignments  of  the  resampled  grid,  each  shifted  by 
0.0 1  cm*1.  The  error  spectra  are  the  difference  between  the  spectra  with  boxcaring  and  w  ithout 
boxcaring  but  resampled  on  the  grid  of  the  boxcared  spectrum.  The  error  spectra  are  shown  both 
with  and  without  deconvolution. 

Considering  the  deconvolved  spectra,  the  maximum  error  is  as  small  as  0.05  percent,  in 
the  case  where  the  line  center  falls  on  a  resampled  grid  point  (solid  line),  and  as  large  as  0.2 
percent  when  the  line  center  falls  halfway  between  the  resampled  grid  points  (dashed  line). 
Without  deconvolution,  the  corresponding  errors  are  0.2  and  0.3  percent  respectively.  Since  the 
alignment  of  spectral  features  relative  to  the  sampling  grid  is  arbitrary,  the  maximum  error  in 
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Figure  4.  An  Example  of  the  Error,  With  and  Without  Deconvolution,  Introduced  by  Boxcaring 
the  Input  Spectrum.  The  input  spectrum  consists  of  a  single  Lorenz  line  of  HWHM  of  0.04 
cm'*  on  a  0,01  cm'*  grid.  The  scanning  function  is  a  sine  of  HWHM  of  0.3  cm'*,  with  r  -  8 
and  M  =  7,  and  the  output  grid  spacing  is  0.07  cm1.  The  four  lines  on  each  error  plot 
correspond  to  different  allignments  of  the  output  grid  relative  to  the  input  grid  (see  text.) 


each  case  must  be  assumed.  Under  this  assumption,  the  deconvolution  process  decreases  the 
boxcaring  error  by  only  one  third.  However,  deconvolution  represents  only  a  small  percent  of 
the  computational  time  so  that  deconvolution  is  still  worth  the  computational  cost. 

The  level  of  error  seen  in  this  figure  is  typical  of  the  error  for  a  value  of  MRATIO  of  8. 
However,  the  actual  error  in  a  particular  situation  will  vary  from  this  example.  For  cases  where 
high  accuracy  (error  less  than  0.5  percent)  is  required,  the  user  is  urged  to  experiment  with 
different  values  of  the  parameters  MRATIO  and  IVX  and  verify  the  accuracy  for  the  particular 
application.  The  default  value  of  MRATIO  adopted  for  FFTSCAN  is  12. 


1 1 


3.  Fast  Fourier  Transform 


The  basic  Fast  Fourier  Transform  package  used  here  is  from  Press5.  Since  the  spectrum 
is  a  real  function,  its  transform  is  Hermitian,  that  is,  the  negative  components  are  the  complex 
conjugate  of  the  positive  components  and  only  the  positive  components  need  be  stored.  We  used 
the  subroutine  REALFT.  FOR,  which  can  calculate  both  the  forward  transform  of  a  real  function 
to  produce  an  Hermitian  function,  and  the  inverse  transform  of  an  Hermitian  function  to  produce 
a  real  function. 

A  FASCODE  spectrum  can  contain  many  thousands  of  points,  more  than  can  be  stored  in 
memory  and  transformed  in  place.  Mark  Esplin,  of  Stewart  Radiance  Lab  has  kindly  provided  a 
disk-based  FFT  routine  which  can  transform  an  array  of  any  size6.  In  the  disk  swapping  FFT. 
the  set  of  input  data  points  is  divided  into  blocks  and  written  to  a  direct  access  file.  Only  two  of 
these  blocks  of  data  reside  in  the  central  memory  of  the  computer  at  a  given  time.  As  the 
transformation  progress,  these  blocks  of  data  are  read,  processed,  and  then  rewritten  to  the  same 
locations.  As  the  routine  proceeds,  the  data  is  First  sorted  into  a  particular  order,  an  in-memory 
FFT  (the  same  REALFT .  FOR  mentioned  previously)  is  then  applied  to  each  block  of  data,  and 
finally  the  data  from  the  various  blocks  are  combined  to  form  the  Fourier  transformation  of  the 
entire  data  set.  The  manner  in  which  the  data  are  sorted  into  the  appropriate  blocks  and  the  way 
the  data  from  the  blocks  are  combined  into  the  Fourier  transformation  of  the  entire  data  set  is  an¬ 
alogous  to  that  of  the  standard  FFT.  In  this  analogy,  blocks  of  data  correspond  to  the  individual 
elements  of  the  standard  FFT  and  an  array  of  data  blocks  on  the  mass  storage  device  corresponds 
to  the  linear  array  of  input  data  points.  In  addition  to  the  blocks  of  data,  there  are  disk  blocks 
generated  that  contain  sine-cosine  information.  The  number  of  these  blocks  is  1/8  the  number  of 
data  blocks.  Both  the  size  of  each  block  and  the  number  of  data  blocks  need  be  a  power  of  2. 

The  parameter  LPTSMX  determines  the  block  size  for  the  FFT.  (See  Section  6  for  a 
discussion  of  this  parameter.)  If  the  spectrum  has  fewer  than  LPTSMX  points,  then  the  FFT  is 
done  in  memory:  if  it  has  more,  the  disk-based  FFT  is  used.  In  either  case,  the  spectrum  is  zero 
filled  as  needed  so  that  the  total  number  of  points  (and  of  blocks  for  the  disk-based  FFT)  is  a 
power  of  2.  The  difference  in  computational  time  for  the  in-memory  versus  disk-based  FFT  for 
the  same  number  of  points  varies  depending  on  the  target  computer.  The  following  table  gives 
one  example  of  a  spectrum  of  131072  =  217  points  scanned  both  ways  on  a  Sun  SparcStation.  In 
this  case,  the  disk-based  case  took  almost  3  times  as  long  as  the  in-memory  case.  On  other 
machines,  the  ratio  is  reported  to  be  smaller. 


LPTSMX  Blocks  Execution  Time  (sec) 

In  Memory 

Disk-Based 

131072  l  21.5 

4096  32  55.5 

^  W.  H.  Press,  B.  P.  Flannery,  S.  A.  Teukolsky,  and  W.  T.  Vetterling,  1987:  Numerical  Recipes 
in  FORTRAN ,  Cambridge  University  Press,  NY. 

6  Mark  Esplin,  Stewart  Radiance  Lab,  Bedford,  MA.,  private  communication. 


4.  Program  Instructions  and  Notes 

4.1  User  Instructions 

FFTSCAN  is  controlled  by  a  single  input  recoid,  similar  to  that  for  the  normal  scanning 
function.  The  details  of  this  record  are  as  follows: 


Field: 

HWHM 

VI 

V2 

JEMIT 

JFN 

MRAT 

DVOUT 

Column: 

1-10 

11-20 

21-30 

31-35 

36-40 

41-45 

46-55 

Field: 

IUNIT 

FIL 

NFIL 

JUNIT 

IVX 

NOFIX 

Column: 

56-60 

61-65 

66-70 

71-75 

76-78 

79-80 

Format (3F10. 3, 315 

,  F10 . 3  , 

,  415,13, 

12) 

HWHM  Half  Width  at  Half  Maximum  of  the  scanning  function,  or  if  JFN  <  0.  the 

maximum  optical  path  difference  of  an  equivalent  interferometer.  If  HWHM  <  0. 
then  exit  FFTSCAN. 

VI  Initial  wavenumber  for  the  scanned  result 

V2  Final  wavenumber  for  the  scanned  result 

JEMIT  =  0:  convolve  with  transmittance 
=  1 :  convolve  with  radiance 

JFN  Selects  the  Scanning  Function  (See  Table  l) 

=  0:  Boxcar.  Halfwidth  is  truncated  to  M  do/2,  where  M  is  an  integer  and  do  is 
the  grid  spacing  of  the  input  spectrum 
=  1:  Triangle 
=  2:  Gauss 
=  3:  Sine2 
=  4:  Sine 
=  5:  Beer 
=  6:  Hamming 
=  7:  Hanning 

If  JFN  <  0,  then  HWHM  is  the  maximum  optical  path  difference  of  an  equivalent 
interferometer,  apodized  to  give  the  scanning  function  given  by  IJFNI. 

MRAT  For  prescanning  with  a  boxcar,  the  ratio  of  HWHM  of  the  scanning  function  to 
the  halfwidth  of  the  boxcar,  default  =12.  If  MRAT  <  0,  no  boxcaring  is 
performed  (see  Notes.) 

DVOUT  Output  grid  spacing  (Not  used,  reserved  for  future  use) 
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IUNIT  Unit  number  of  the  tile  containing  the  spectrum  to  be  scanned,  default  =  i  1  (see 
Notes) 

IFIL  Sequential  number  of  the  first  FASCODE  file  on  IUNIT  to  be  scanned 

NFIL  Number  of  FASCODE  files  on  IUNIT  to  be  scanned, 
beginning  with  IFIL 

JUNIT  Unit  number  of  the  file  containing  the  output  spectrum,  default  =  12  {see  Notes  ) 

IVX  =  -1:  Scanning  function  is  calculated  as  the  FFT  of  the  Apodization  function 
=  0:  Program  decides  how  to  calculate  the  scanning  function,  using  CR  from 
Table  1. 

=  1:  Scanning  function  is  calculated  analytically 

NOFIX  For  prescanning  with  a  boxcar:  if  non-zero,  then  do  not  deconvolve  the  scanned 
spectrum  with  the  boxcar 


4.2  Notes 

The  program  expands  the  spectral  interval  V]_  to  V2  to  V]_  -  Cr(JFN)*HWHM  to  V2  + 
Cr(JFN)*HWHM  before  smoothing  (see  Table  1.)  This  expansion  ensures  that  edge  effect  do 
not  contaminate  the  endpoints  of  the  scanned  spectrum.  If  there  is  insufficient  data  for  this 
expansion,  then  the  maximum  expansion  possible  is  performed  and  an  informative  message  is 
written.  The  output  spectrum  extends  from  V^-2d\)  to  V2+2do,  where  do  is  the  output  spacing. 
The  extra  points  at  either  end  allow  for  four-point  interpolation  of  the  spectrum  at  the  endpoints. 

The  default  value  of  MRAT  should  reduce  the  boxcaring  error  sufficiently  for  most 
applications  (better  than  0.2  percent).  For  greater  accuracy,  it  is  necessary  to  increase  MRAT, 
or  set  it  to  -1  to  turn  off  boxcaring. 

The  spectral  input  and  output  files  are  on  units  IUNIT  and  JUNIT,  which  default  to  12 
and  1 1  respectively.  For  input,  if  no  file  is  open  on  unit  =  IUNIT  then  the  program  looks  for  a 
file  named  TAPExx,  where  xx  =  IUNIT  (e.g.  TAPE12.)  If  a  file  by  this  name  does  not  exist, 
then  an  error  results  and  then  the  program  stops.  For  output,  if  a  file  is  open  on  unit  =  JUNIT  . 
then  that  file  is  rewound  and  used.  If  not,  then  the  program  looks  for  a  file  by  named  TAPExx . 
If  that  file  does  exist,  then  an  error  results  and  the  program  stops.  If  it  does  not  exist,  then  a  new 
file  by  that  name  is  opened  for  output.  If  IUNIT  or  JUNIT  are  negative,  then  the  program 
reads  in  the  spectral  file  names  (60  characters  maximum,  including  the  path)  on  the  next  record. 
If  the  named  input  file  does  not  exist  or  if  the  named  output  file  does  exist,  then  an  error  results 
and  the  program  stops.  Otherwise,  they  are  opened  on  the  first  free  unit  numbers  between  61 
and  99. 
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5.  Examples 

Figure  5  shows  an  example  of  a  FASCODE  spectrum  smoothed  using  FFTSCAN.  The 
calculated  spectrum  models  the  upward  radiance  at  72  km  for  the  US  Standard  Atmosphere. 
The  monochromatic  calculation  extended  from  600  to  800  cm'1  with  a  do  of  0.000953.  In 
Figure  5a,  the  monochromatic  FASCODE  spectrum  was  smoothed  with  a  sine  scanning  function 
of  HWHM  =  0.2196  cm'1,  corresponding  to  the  HIS  resolution.  Figure  5b  shows  the  error  in  the 
scanned  spectrum  from  using  boxcaring  with  deconvolution  (r  =  12).  This  value  of  r  (or  MRAT 
in  the  user  instructions)  gives  a  value  of  M  of  38,  resulting  in  a  38  fold  reduction  in  the  number 
of  points  in  the  spectrum.  The  maximum  error  of  about  lxl0‘8  is  about  0.5  percent  of  the 
typical  spectral  excursion  of  2xl0‘6  or  about  0.2  percent  of  the  maximum  spectral  value  of 
8xlO'8.  For  reference.  Figure  5c  shows  the  error  using  the  standard  FASCODE  3  convolution 
with  a  bound  of  80  halfwidths.  The  reference  spectrum  for  calculating  the  errors  in  Figures  5b 
and  c  is  the  FFTSCAN  calculation  without  boxcaring,  interpolated  to  the  proper  grid. 

Table  2  compares  the  computational  time  and  maximum  error  for  the  calculations  shown 
in  Figure  5.  The  calculations  were  performed  on  an  Sun  SparcStation  and  the  scanning 
functions  were  applied  from  625  cm-1  to  775  cm'1.  The  results  for  scanning  with  the  sine-  and 
the  triangular  scanning  functions  are  also  shown.  Note  that  the  boxcaring  errors  for  the  other 
functions  are  four  times  less  than  that  for  the  sine  function. 

These  results  show  that  for  the  sine  function,  FFTSCAN  with  boxcaring  provides  a  three¬ 
fold  increase  in  speed  and  a  better  than  two  fold  increase  in  accuracy  over  FASCODE.  For  the 
sine2  and  the  triangle,  the  execution  time  for  the  two  programs  is  about  equal,  but  FFTSCAN  is 
again  twice  as  accurate.  The  reference  spectrum  for  calculating  the  errors  is  again  the 
FFTSCAN  calculation  without  boxcaring,  interpolated  to  the  proper  grid. 


15 


650  660  670  680  690  700 


Figure  5.  Example  of  a  Smoothed  Spectrum  and  the  Associated  Error:  a.  FASCODE  Calculated 
Spectrum  Smoothed  by  FFTSCAN,  b.  The  Error  In  the  Smoothed  Spectrum  from  Using 
FFTSCAN  With  Boxcaring,  and  c.  Error  From  Using  the  FASCODE3  Scanning  Function.  The 
scanning  function  is  a  sine  with  a  HWHH  of  0.2196  cm'1,  corresponding  to  the  HIS  instrument 
in  the  unapodized  mode. 
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Table  2.  Comparison  of  Computational  Time  and  Accuracy,  FFTSCAN  versus  FASCODE, 
For  the  Sine,  Sine2,  and  Triangle  Scanning  Functions. 


Computational  Time 
(sec) 

Maximum  Error 
(Percent  Radiance) 

Scanning  Option 

Sine 

Sine2 

Triangle  ! 

Sine 

Sine2 

Triangle 

FFTSCAN,  No  Boxcar 

117 

117 

117 

(NA) 

(NA) 

(NA) 

FFTSCAN,  With  Boxcar 

5.8 

3.7 

3.3 

0.2 

0.05 

0.05 

FASCODE3  Convolution 

3.9 

3.0 

0.5 

0.1 

0.1 

Notes: 

a.  The  computational  times  refer  to  a  Sun  SparcStation.  Times  are  approximate  and  both 
the  absolute  and  the  relative  times  will  vary  depending  on  the  case. 

b.  The  spectral  extent  of  the  smoothed  spectrum  was  from  650  to  775  cm-1 

c.  The  monochromatic  do  was  0.000953  env1 

d.  The  number  of  points  in  the  scanned  function  was  167,958  (noboxcaring)  and  4419  (with 
boxcaring),  a  reduction  of  a  factor  of  38. 
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6.  Implementation  Notes 

FFTSCAN  is  written  in  ANSI  Standard  Fortran  77  and  is  designed  to  be  highly  portable. 
It  was  developed  on  Sun  workstation  under  Unix,  and  early  versions  have  been  ported  to  a 
VAX,  a  Cyber  computer  under  NOS/VE  and  an  IBM  PC.  FFTSCAN  is  designed  either  to  be  run 
as  an  independent  program  or  to  be  included  as  a  module  of  FASCODE. 

The  program  uses  three  include  files:  f  f tparm.  inc ,  parmcomm . inc ,  and 
scancomm.  inc.  These  files  carry  parameters,  common  blocks,  and  other  declarations  used 
throughout  the  program.  The  file  scancomm.  inc  includes  the  DOUBLE  PRECISION 
statements: 

Implicit  Double  Precision  (V) 

Double  Precision  XID, SECANT, HMOLID, XALTZ, YID 

which  may  or  may  not  have  to  be  disabled,  depending  on  how  it  is  set  in  FASCODE.  Typically, 
it  is  enabled  on  a  32  bit  machine  and  disabled  on  a  64  bit  machine.  There  is  no  standard  syntax 
for  the  include  statement  and  these  statements  may  have  to  be  changed  to  suit  a  particular 
system.  For  versions  of  Fortran  which  do  not  support  include  files,  the  program  may  be 
distributed  with  these  files  already  included. 

There  are  a  few  hardware  dependencies  related  to  the  disk-ba*w»d  FFT  which  must  be 
considered.  If  the  size  of  the  spectrum  is  greater  than  memory'  set  aside  for  the  in-memory  FFT. 
then  the  program  uses  the  disk-based  FFT.  The  disk-based  FFT  divides  the  input  data  points 
into  blocks  and  writes  these  blocks  to  disk  as  direct  access  records.  An  in-memory  FFT  is  per¬ 
formed  on  each  record.  The  size  of  each  record  must  be  a  power  of  2.  The  parameter  LPTSMX 
sets  the  maximum  size  of  an  in-me  ory  FFT  and  the  size  of  the  direct  access  records.  This  vari¬ 
able  should  be  set  as  large  as  possible  for  the  particular  computer  since  the  in-memory  FFT  is 
more  efficient.  The  pitfall  that  must  be  avoided  is  setting  LPTSMX  too  large  in  which  case  the 
virtual  memory  system  will  page  the  data.  Since  the  points  processed  by  the  in-memory  FFT 
come  from  widely  scattered  locations,  the  number  of  calls  to  the  disk  will  be  extremely  large 
(thrashing)  and  the  performance  will  be  very  poor.  On  the  other  hand,  the  minimum  size  of  an 
FFT  is  also  LPTSMX  (this  is  a  design  error  which  will  be  corrected  in  the  next  version.)  If 
LPTSMX  is  large  but  the  region  to  be  scanned  is  small,  then  the  calculation  will  take  unnecessar¬ 
ily  long.  Therefore,  LPTSMX  should  be  set  to  somewhere  between  the  smallest  typical  spectral 
size  and  the  largest  value  possible  without  thrashing.  The  parameter  LPTSMX  is  set  in  the 
include  file  ff  tparm.  inc.  The  user  may  have  to  adjust  LPTSMX  from  the  default  value 
(8192  =  213)  to  be  optimum  for  the  target  computer. 

To  determine  whether  the  program  is  using  the  in-memory  or  the  disk-based  FFT,  note 
the  following  line  in  the  output: 

FFT:  Total  number  of  points  =  xxxxx  and  blocks  =  vy 

If  the  number  of  points  yy  is  2  or  greater,  then  the  disk-based  FFT  is  being  used. 

The  parameter  IBLKSZ,  declared  in  the  same  statement,  is  the  block  size  of  the  direct 
access  records  and  is  used  in  OPEN  statements.  Depending  on  the  computer,  it  may  be  in  words 
(e.g.,  VAX  and  CDC  Cyber)  in  which  case  IBLKSZ  =  LPTSMX,  or  in  bytes  (Microsoft  Fortran 
for  the  PC,  Sun,  and  Alliant)  in  which  case  IBLKSZ  =  4  *LPTSMX.  On  a  CRAY,  which  has  a 
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word  size  of  64  bits  (B  bytes),  use  IBLKSZ  =  8*LPTSMX.  In  VAX  Fortran,  the  maximum 
direct  access  block  size  is  4095  words  so  that  the  maximum  allowed  value  of  IBLKSZ  is  2048. 
Since  this  is  much  smaller  than  the  maximum  physical  memory,  the  user  may  want  to  use  a  large 
value  of  LPTSMX  and  disable  the  disk-based  FFT  capability. 

On  CDC  Cyber  machines,  the  subroutines  BUFIN  and  BUFOUT  may  have  to  be  changed 
to  match  the  those  used  in  FASCODE. 

To  modify  FFTSCA N  to  include  other  scanning  functions,  the  subroutines  FFTSCN.F 
and  SCNFNT.F  need  to  be  modified.  In  FFTSCN.F,  modify  the  variables  JFNMAX. 
ANAMSES,  C,  CRATIO,  and  CLIMIT  as  appropriate  (the  definitions  of  these  variables  are 
given  in  program  comments.)  In  SCNFNT.F,  use  an  existing  scanning  function  as  a  example, 
and  add  the  equations  defining  the  new  function,  both  in  the  frequency  and  in  the  space  domain. 


19 


References 


Clough,  S.,  A.,  F.  X.  Kneizys,  L.  S.  Rothman,  and  W.  O.  Gallery,  1981, Atmospheric  spectral 
transmittance  and  radiance:  FASCOD1B,  Proc.  of  Soc.  Photo.  Opt.  Instrm.  Ena.,  277, 
152-166. 

Clough,  S.  A.,  F.  X.  Kneizys,  G.  P.  Anderson.  E.  P.  Shettle.  J.  H.  Chetwynd.  L.  W.  Abreu.  and 
L.  A.  Hall,  1989,  FASCOD3:  Spectral  Simulation,  in  IRS  '88:  Current  Problems  in 
Atmospheric  Radiation.  J.  Lenobie  and  J.  F.  Gelyn  (Eds.),  A.  Deepak  Pub. 

Bell,  J.  R.,  1972,  Introductory  Fourier  Transform  Spectroscopy.  Academic  Press,  New  York 

Esplin,  Mark,  Stewart  Radiance  Lab,  Bedford,  MA.,  private  communication. 

Press,  W.  H.,  B.  P.  Flannery,  S.  A.  Teukolsky,  and  W,  T.  Vetterling,  1987:  Numerical  Recipes  in 
FORTRAN,  Cambridge  University  Press,  NY. 

Revercomb,  H.  E..  H.  Buijs.  H.  B.  Howell,  D.  D.  LaPorte.  W.  L.  Smith,  and  L.  A.  Sromovsky. 
1988,  "Radiometric  Calibration  of  IR  Fourier  Transform  Spectrometers:  Solution  to  a 
Problem  with  the  High-Resolution  Interferometer  Sounder,"  Appl.  Opt.  27,  p32 10 


20 


